Journal  of  Power  Sou 


;  215  (2012)  135-144 


Contents  lists  available  at  SciVerse  ScienceDirect 

Journal  of  Power  Sources 


ELSEVIER 


epage:  www.e 


m/locate/j  powsou  r 


A  mathematical  model  for  predicting  the  life  of  polymer  electrolyte  fuel  cell 
membranes  subjected  to  hydration  cycling 

S.F.  Burlatsky3,  M.  Gummalla3  *,  J.  O’Neill d,  V.V.  Atrazhevb  c,  A.N.  Varyukhinc,  D.V.  Dmitriev  b  c, 
N.S.  Erikhman  b,c 

a  United  Technologies  Research  Center,  411  Silver  Lane,  East  Hartford,  CT  06108,  USA 
b  Russian  Academy  of  Science,  Institute  of  Biochemical  Physics,  Kosygin  str.  4,  Moscow  119334,  Russia 
c  Science  for  Technology  LLC,  Leninskiy  pr-t  95,  Moscow  119313,  Russia 
d  UTC  Power,  195  Governor’s  Highway,  South  Windsor,  CT  06074,  USA 


HIGHLIGHTS 


►  A  novel  modeling  framework  to  predict  the  PEMFC  membrane  lifetime. 

►  Membrane  lifetime  as  a  function  of  RH  cycling  amplitude  and  membrane  mechanical  properties. 

►  Very  good  agreement  with  the  experiments. 

►  Membrane  failure  in  the  fuel  cell  is  caused  by  damage  accumulation  under  cyclic  stress  in  the  membrane. 

►  Tensile  stress  causes  a  crack  formation  at  RH  cycled  side  of  the  membrane. 
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Under  typical  Polymer  Electrolyte  Membrane  Fuel  Cell  (PEMFC)  fuel  cell  operating  conditions,  part  of  the 
membrane  electrode  assembly  is  subjected  to  humidity  cycling  due  to  variation  of  inlet  gas  RH  and/or 
flow  rate.  Cyclic  membrane  hydration/dehydration  would  cause  cyclic  swelling/shrinking  of  the 
unconstrained  membrane.  In  a  constrained  membrane,  it  causes  cyclic  stress  resulting  in  mechanical 
failure  in  the  area  adjacent  to  the  gas  inlet.  A  mathematical  modeling  framework  for  prediction  of  the 
lifetime  of  a  PEMFC  membrane  subjected  to  hydration  cycling  is  developed  in  this  paper.  The  model 
predicts  membrane  lifetime  as  a  function  of  RH  cycling  amplitude  and  membrane  mechanical  properties. 
The  modeling  framework  consists  of  three  model  components:  a  fuel  cell  RH  distribution  model, 
a  hydration/dehydration  induced  stress  model  that  predicts  stress  distribution  in  the  membrane,  and 
a  damage  accrual  model  that  predicts  membrane  lifetime.  Short  descriptions  of  the  model  components 
along  with  overall  framework  are  presented  in  the  paper.  The  model  was  used  for  lifetime  prediction  of 
a  GORE-SELECT  membrane. 

©  2012  Published  by  Elsevier  B.V. 


1.  Introduction 

Mechanical  degradation  of  PEM  membrane  can  limit  stack  life¬ 
time.  Operation  of  a  fuel  cell  under  realistic  load  cycles  results  in 
both  chemical  and  mechanical  degradation  of  the  polymer  elec¬ 
trolyte  membrane.  Degradation  of  the  membrane  causes  opening 
up  of  pinholes  or  crazing  of  polymer  [1—3]  increasing  gas  crossover 
and  subsequently  resulting  in  catastrophic  failures  of  the  fuel  cell 
stack.  The  review  article  by  Borup  et  al.,  captures  the  degradation 
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mechanisms  in  the  PEM  fuel  cells,  which  includes  both  membrane 
and  electrode  degradation.  The  synopsis  of  the  membrane  degra¬ 
dation  mechanisms  is  provided  with  radical  attack  of  the  polymer 
chains  being  the  primary  cause  for  chemical  degradation.  Due  to 
constant  compressive  forces  on  the  polymer  associated  with  the 
cell  build  the  polymers  undergo  physical  degradation,  resulting  in 
polymer  creep.  It  is  suggested  that  the  physical  degradation  modes 
are  relatively  less  studied.  However,  post  test  analysis  of  the 
membranes  suggest  localized  failures  in  the  membranes  at  the  high 
stress  points,  such  as  the  edge  of  the  ribs  [4], 

Understanding  and  modeling  the  mechanical/physical  degra¬ 
dation  mechanism  and  kinetics  enables  prediction  of  membrane 
lifetime  as  a  function  of  PEM  operational  conditions  and 
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optimization  of  membrane  structure  through  the  choice  of  rein¬ 
forcement  [5—8],  the  membrane  processing  methods  and  the 
operating  conditions.  The  physics  based  model  of  membrane 
mechanical  degradation  could  guide  the  synthesis  of  new 
membranes  with  tuned  mechanical  properties  and  enhanced  life  in 
a  fuel  cell. 

Hydration  cycling  is  a  primary  cause  of  the  mechanical  degra¬ 
dation  of  a  geometrically  constrained  polymer,  which  exhibits 
dimensional  changes  with  varied  water  content.  A  simple  way  to 
impose  mechanical  damage  to  a  geometrically  constrained  PEM 
membrane  is  to  subject  it  to  humidity  cycling.  At  high  RH,  the 
membrane  absorbs  water,  and  at  low  RH,  the  membrane  desorbs 
water.  Such  RH  cycling  would  result  in  swelling  and  shrinking  of  an 
unconstrained  polymer.  RH  cycling  of  a  constrained  polymer  cau¬ 
ses  cyclic  stress.  In  PEM,  such  geometrical  constraints  are  imposed 
by  bipolar  plate  ribs  through  the  gas  diffusion  layers  (GDLs), 
catalyst  layers  adjacent  to  the  membrane  and  at  the  seals.  More¬ 
over,  we  hypothesize  that  internal  stress  in  the  membrane  can  be 
caused  by  the  difference  in  gas  RH  at  the  anode  and  cathode 
membrane  surfaces  that  routinely  occurs  at  fuel  cell  operational 
conditions.  Cyclic  stress  in  the  membrane  causes  irreversible 
elongation  of  the  membrane  [3]  and  subsequent  formation  of 
crazes  and  cracks  that  causes  gas  crossover  through  the  membrane 
and  stack  failure. 

The  goal  of  this  work  is  to  develop  a  model  that  predicts 
membrane  lifetime  as  a  function  of  fuel  cell  design  and  operating 
conditions,  and  also  of  membrane  transport  and  mechanical 
properties.  Three  components  are  needed  to  model  the  mechanical 
degradation  process  under  RH  cycling.  These  components  are: 
a  fuel  cell  RH  distribution  model,  a  hydration/dehydration  induced 
stress  model,  and  a  damage  accrual  model,  see  Fig.  1.  The  fuel  cell 
RH  distribution  model  calculates  RH  distribution  in  fuel  cell  gas 
channels  as  a  function  of  operating  conditions  and  time.  This 
distribution  depends  on  fuel  cell  design.  The  stress  model  predicts 
the  stress  profile  in  the  membrane  for  a  given  time-dependent 
profile  of  RH  at  the  membrane/electrode  interfaces.  The  damage 
accrual  model  predicts  the  membrane  lifetime  for  a  given  stress 
profile  in  the  membrane.  The  latter  two  model  components  use  the 
extended  Eyring  model  of  polymer  viscoelastic  deformation.  The 
damage  accrual  model  component  predicts  the  membrane  irre¬ 
versible  elongation  as  a  function  of  applied  stress  and  time.  RH  and 
temperature  distribution  in  gas  channels  and  membrane  hydration 
at  steady-state  conditions  as  a  function  of  fuel  cell  operating 
conditions  were  published  and  discussed  in  [9—15],  Hence,  they  are 
not  described  further  in  this  work,  and  the  focus  is  on  other  model 
components.  The  main  features  and  limitations  of  the  three  novel 
modeling  components  are  briefly  discussed  below. 

f.f.  Viscoelastic  polymer  deformation 

A  large  number  of  models  of  polymer  deformation  are  available 
in  the  literature.  The  linear  theory  of  viscoelasticity  was  introduced 
by  Boltzmann  [16]  many  decades  ago,  and  it  provides  the  basis  for 
all  well-known  constitutive  models  of  linear  viscoelasticity 


(Maxwell  model,  Kelvin— Voigt  model,  Standard  linear  solid  model, 
and  their  generalizations).  These  models  utilize  the  analogy  with 
mechanical  systems  consisting  of  elastic  springs  and  viscous 
dampers.  Models  of  linear  viscoelasticity  often  fail  when  either 
high  deformation  (>10%)  or  long-term  behavior  of  polymers  is 
studied.  Therefore,  to  predict  membrane  lifetime,  an  approach  is 
required  which  accounts  for  non-linear  effects  of  polymer 
viscoelasticity. 

There  are  several  approaches  for  non-linear  viscoelasticity 
modeling.  The  empirical  approach  uses  correlations  between  time- 
dependent  stress  and  strain  [17—19],  This  approach  relies  on 
numerous  fitting  parameters  that  are  specific  for  a  given  polymer. 
The  alternative  approach  is  a  semi-empirical  [20,21]  or  purely 
mathematical  [22]  generalization  of  Boltzmann’s  linear  theory.  The 
non-linear,  time-dependent  constitutive  model  for  prediction  of 
the  hygro-thermomechanical  behavior  of  Nafion  was  proposed  in 
[23]  and  fitted  to  experimental  data.  Neither  approach  accounts  for 
the  microscopic  physics  of  polymer  deformation.  In  contrast  to  the 
empirical  and  semi-empirical  models,  the  molecular  theory  of  non¬ 
linear  viscoelasticity  proposed  by  Eyring  et  al.  [24]  is  based  on 
a  physical  concept  of  polymer  dynamics.  This  model  assumes  that 
polymer  deformation  occurs  as  a  motion  of  polymer  chain 
segments  that  overcome  potential  barriers  at  the  entanglement 
points.  This  model  predicts  the  (non-linear)  dependence  for  poly¬ 
mer  elongation  on  applied  constant  stress.  In  this  paper,  the  Eyring 
concept  is  extended  to  predict  stress  relaxation  in  a  constrained 
polymer. 

1.2.  Hydration/dehydration  stress  model 

The  hydration/dehydration  stress  model  calculates  water 
distribution  in  the  membrane  as  a  function  of  time  and  local  stress 
in  the  membrane  caused  by  changes  in  water  content.  The  equi¬ 
librium  water  content  of  Nafion  membranes  as  a  function  of  gas  RH 
was  experimentally  studied  in  [25].  The  dynamics  of  water  sorp¬ 
tion/desorption  by  Nafion  membranes  and  water  transport  through 
membranes  were  studied  in  [26—28],  Nafion  demonstrates  unusual 
water  transport  properties.  For  example,  its  water  sorption  time  is 
an  order  of  magnitude  larger  than  its  water  desorption  time  [26], 
Several  mathematical  models  of  water  transport  in  Nafion 
membranes  are  available  in  the  literature  [29-32],  Currently, 
models  of  ionomer  water  transport  are  based  on  the  diffusion 
equation,  and  modeling  efforts  are  focused  on  calculation  of  the 
water  diffusion  coefficient  as  a  function  of  membrane  water 
content.  The  ionomer/gas  interfacial  barrier  was  hypothesized  in 
[32],  The  model  developed  in  [32]  explains  the  peculiar  kinetics  of 
water  sorption/desorption  by  the  kinetics  of  water  evaporation/ 
condensation  at  the  ionomer/gas  boundary.  In  the  current  paper, 
we  assume  thermodynamic  equilibrium  of  water  in  the  membrane 
with  vapor  at  the  membrane  interface,  taking  advantage  of  the  low 
ratio  of  water  sorption/desorption  time  constant  to  the  humidity 
cycling  period.  Water  transport  in  the  membrane  bulk  is  treated  as 
diffusion,  with  the  water  diffusion  coefficient  dependent  on  the 
membrane  water  content. 


Membrane  Mechanical  Degradation  Model: 


Fig.  1.  Schematic  diagram  showing  the  components  of  the  model  framework  to  predict  mechanical  life  of  the  membrane  undergoing  hydration  cycles. 
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The  hydration/dehydration  stress  model  assumes  that  the  local 
stress  in  the  membrane  is  induced  by  competition  between  the 
membrane’s  tendency  to  swell  or  shrink  in  response  to  a  change  in 
hydration  level  and  the  geometrical  constraints  that  prevent  the 
membrane  from  swelling  or  shrinking.  Unconstrained  PFSA 
membranes  are  known  to  absorb  water  and  undergo  dimensional 
changes  [1],  The  dimensional  changes  can  be  controlled  by  the 
choice  of  reinforcement  [6]  and  the  membrane  processing 
methods.  Reinforced  membrane  demonstrates  approximately  4x 
lower  dimensional  change  than  that  of  non-reinforced  membrane 
under  the  same  conditions.  Tang  et  al.  experimentally  studied 
mechanical  properties  and  dimensional  change  of  Nafion  111 
membrane  [1]  and  reinforced  Gore  membrane  [6]  as  a  function  of 
temperature  and  RH.  According  to  [1,6]  the  membrane  dimensional 
change  is  proportional  to  the  change  of  RH  in  ambient  air.  Though, 
the  dependence  of  water  content  in  the  membrane  on  RH  [33]  is 
a  non-linear  function,  we  linearized  it  in  the  water  content  interval 
from  \  =  3  to  \  =  12  to  speed  up  the  calculations.  In  fact,  we  carried 
out  some  numerical  calculations  with  more  accurate  approxima¬ 
tion  X(RH)  and  found  only  minor  effect  on  the  results.  In  the 
hydration/dehydration  stress  model,  we  utilize  a  new  non-linear 
equation  for  polymer  stress/elongation  that  is  predicted  by  the 
extended  Eyring  model.  Taking  advantage  of  this  stress/elongation 
equation,  along  with  the  linear  dependence  of  the  membrane 
swelling  on  water  content,  we  calculate  the  stress  distribution 
induced  by  cyclic  hydration/dehydration.  The  model  parameters 
are  obtained  from  experimental  data  presented  in  [6,33]. 

1.3.  Damage  accrual  model 

The  damage  accrual  model  is  based  on  experimental  data  for 
membrane  failure  under  cyclic  stress  and  on  our  concept  of  poly¬ 
mer  plastic  deformation.  Tensile  stress  caused  by  membrane 
dehydration  in  a  constrained  Nafion  111  membrane  was  studied 
experimentally  by  Tang  et  al.  [3],  Dependence  of  stress  on  dehy¬ 
dration  level  was  measured  and  the  typical  stress  value  was 
approximately  1  MPa  for  a  membrane  that  was  dehydrated  from 
100%  RH  to  60%  RH.  Such  stress  is  approximately  an  order  of 
magnitude  lower  than  the  tensile  strength  of  a  Nafion  membrane. 
Tang  et  al.  [3]  have  experimentally  shown  that  the  amplitude  of  the 
cyclic  stress  that  caused  substantial  permanent  elongation  of 
Nafion  111  is  1/10  of  the  Nafion  111  tensile  strength.  Nafion 
membrane  creep  tests  under  constant  load  are  reported  in  the  work 
of  Majsztrik  et  al.  [34],  In  [34],  the  Nafion  creep  dynamics  was 
experimentally  studied  as  a  function  of  the  temperature  and  the 
membrane  hydration  level.  The  applied  tensile  stress  was  approx¬ 
imately  equal  to  1.55  MPa.  That  stress  is  also  by  approximately  one 
order  of  magnitude  lower  than  the  tensile  strength  of  Nafion.  A 
strong  dependence  of  Nafion  creep  rate  on  the  hydration  level  was 
observed  in  [34],  At  8%RH  and  >60  °C,  low  water  content  hardens 
the  Nafion  membrane  and  decreases  the  rate  of  creep.  The  creep 
rate  of  dry  Nafion  dramatically  increases  at  ~90  °C,  which  corre¬ 
sponds  to  the  temperature  of  Nafion  a-relaxation  detected  in  DMA 
tests  in  [35,36],  At  the  same  temperature,  the  creep  rate  of  wet 
Nafion  (RH  >  8%)  increases  much  more  slowly  than  that  of  dry 
Nafion.  Authors  explain  observed  Nafion  stiffening  with  RH 
increase  by  the  increased  electrostatic  interaction  in  ionic  clusters 
composed  of  SO3H  polar  groups  in  the  presence  of  water  molecules. 
To  develop  the  microscopic  model  that  explains  experimental  data 
on  Nafion  creep  dynamics,  fundamental  insights  into  the  mecha¬ 
nism  for  polymer  damage  accrual  are  needed. 

Kusoglu  [7]  showed  through  simulations  that  a  hydro-thermal 
loading  under  fuel  cell  operating  conditions  results  in  compres¬ 
sive  stress  in  the  polymer.  According  to  [7],  the  stress  level  exceeds 
the  yield  strength,  causing  permanent  damage  to  the  polymer 


membrane.  The  compressive  stress  causes  extrusion  of  membrane 
materials  from  the  compressed  areas  such  as  the  areas  under  the 
seals.  Apparently,  the  membrane  extrusion  does  not  cause  the  loss 
of  membrane  integrity.  We  assume  that  the  tensile  stress  in 
membrane  is  much  more  damaging  for  membrane  integrity, 
because  it  causes  formation  and  propagation  of  cracks  and  crazes. 
In  the  present  work,  we  focus  on  generation  of  through-plane 
cracks  in  the  membrane  that  result  in  gas  crossover  through  the 
membrane  and  dramatic  performance  loss. 

In  summary,  membrane  hydration  and  dehydration  changes  the 
local  membrane  water  content,  resulting  in  a  stress  cycle  in  the 
polymer  membrane.  A  correlation  between  hydration  level  and 
stress  cycling  is  needed.  Additionally,  to  predict  life  of  the 
membrane,  a  correlation  between  membrane  macroscopic 
mechanical  properties,  stress  cycling,  and  membrane  damage 
accrual  needs  to  be  developed.  The  overall  model  that  predicts 
membrane  lifetime  as  a  function  of  membrane  mechanical  prop¬ 
erties  and  fuel  cell  operating  conditions  could  guide  the  synthesis 
of  new  membranes  with  tuned  mechanical  properties  and 
enhanced  life  in  a  fuel  cell.  Development  of  such  a  modeling 
framework  to  predict  the  life  of  a  membrane  subjected  to  hydration 
cycles  is  the  focus  of  this  work. 

The  rest  of  the  paper  is  organized  as  follows.  The  extended 
Eyring  model  of  polymer  viscoelastic  deformation  is  presented  in 
Section  2.  This  model  is  used  as  a  basis  of  the  membrane  stress  and 
damage  accrual  models  rather  than  as  a  standalone  framework 
component.  A  model  that  predicts  stress  in  the  membrane  during 
RH  cycling  in  gas  channels  for  given  membrane  mechanical  prop¬ 
erties  is  proposed  in  Section  3.  Next,  a  mathematical  model  iden¬ 
tifying  the  functional  form  that  relates  the  stress  cycle  to 
irreversible  damage  of  the  polymer  is  discussed  in  Section  4.  The 
accrual  of  such  irreversible  damage  results  in  subsequent  polymer 
failure,  determining  the  membrane  life  subjected  to  hydration 
cycling.  Input  parameters  of  the  model  are  summarized  in  Section 
5.  Modeling  results  for  a  GORE-SELECT  membrane  and  discussion  of 
the  results  are  presented  in  Section  6.  Conclusions  are  presented  in 
Section  7. 

2.  Extended  Eyring  model  for  polymer  deformation 

The  time-dependent  response  of  the  stress  to  applied  strain  that 
takes  into  account  viscoelastic  relaxation  in  the  polymer  material  is 
needed  for  two  components  of  the  degradation  model.  The  model 
of  polymer  viscoelastic  deformation  utilized  this  paper  is  based  on 
the  assumption  that  polymer  deformation  occurs  through  the 
transport  of  polymer  chains  through  entanglements.  The  early 
Eyring  model  [24]  is  based  on  a  similar  concept.  The  Eyring  model 
calculates  the  non-linear  dependence  of  polymer  deformation  on 
applied  constant  stress,  and  is  relevant  to  creep  experiments.  In 
a  creep  experiment,  a  sample  is  subjected  to  constant  stress  lower 
than  the  yield  stress,  while  the  elongation  is  monitored.  In  this 
section,  we  briefly  describe  a  new,  extended  Eyring  model  that 
calculates  stress  relaxation  in  a  stretched  constrained  polymer 
under  arbitrary  time-dependent  controlled  deformation. 

The  concept  used  in  the  current  model  is  summarized  below. 
Following  [37],  we  model  a  polymer  as  a  set  of  entangled  chains. 
One  chain  wraps  around  another  chain  and  turns  around  at  the 
entanglement  point.  In  the  current  model,  stress  is  transmitted 
from  one  chain  to  another  through  the  entanglements,  which 
secure  the  mechanical  integrity  of  the  polymer.  A  Chain  Segment 
(CS)  is  defined  as  a  fraction  of  the  chain  confined  between  two 
subsequent  entanglements.  The  macroscopic  deformation  of  the 
polymer  is  calculated  from  the  microscopic  elongation  of  CSs 
driven  by  the  changing  distance  between  entanglements.  The 
macroscopic  stress  is  calculated  from  the  microscopic  tension 
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acting  along  the  CSs.  This  concept  to  a  large  extent  is  similar  to  the 
concept  of  interconnected  elastic  springs  in  elasticity  theory.  The  CS 
is  an  analog  of  the  spring;  the  entanglement  is  an  analog  of  the 
point  where  the  springs  are  attached  to  each  other.  The  elastic 
deformation  in  the  current  polymer  model  is  achieved  through  the 
elastic  deformation  of  individual  CSs,  which  is  an  analog  of  the 
elastic  deformation  of  the  springs  in  elasticity  theory.  The  funda¬ 
mental  difference  comes  from  the  fact  that  the  chains  can  slip 
through  entanglements,  changing  the  equilibrium  length  of  CSs. 
We  assume  that  slippage  of  the  chain  through  entanglements 
causes  irreversible  elongation  of  the  polymer. 

To  calculate  the  slip  rate  as  a  function  of  tensile  force  we 
consider  the  dynamics  of  CS  at  the  entanglement,  as  indicated  in 
Fig.  2a.  In  the  absence  of  the  tensile  force,  F,  the  monomer  adjacent 
to  the  entanglement  is  located  in  the  local  energy  minimum 
between  two  symmetric  energy  barriers  Uo  (Fig.  2b).  The  chain  slips 
through  the  entanglement  when  the  monomer  overcomes  one  of 
the  energy  barriers. 

In  a  stressed  polymer,  a  tensile  force,  F,  is  induced  by  stress  and 
acts  on  the  chain  (Fig.  2a),  disturbing  the  symmetry  of  the  mono¬ 
mer’s  potential  energy  well:  the  right  barrier  becomes  lower  by  Fa 
and  the  left  barrier  becomes  higher  by  Fa,  where  a  is  the  monomer 
length.  This  results  in  preferable  displacement  of  the  monomer  to 
the  right,  i.e.  along  the  force  direction.  Here  we  assume  that  the 
monomer  motion  is  thermally  induced  and  assisted  by  the  force  F, 
Fa  «  Uq.  The  rate  V4"  at  which  the  monomer  jumps  to  the  right  is 


Here  w0~1012  s-1  is  atomic  frequency,  k  =  1.3810-23  J  K-1  is 
Boltzmann  constant  and  T  is  absolute  temperature. 

The  rate  of  monomer  displacement  to  the  left  V~  is 


(2) 

Therefore,  the  average  rate  of  polymer  transport  through  the 
entanglement  is 

V  =  awsinh  (3) 

Here  w  is  the  frequency  of  thermally  activated  jumps  of 
a  monomer  in  the  local  minimum  in  the  entangled  state: 

M  =  w0exp(  ^  (4) 


Here  a2«/kT  is  a  mobility  of  the  chain  in  the  entanglement.  In  the 
opposite  limit,  Fa/kT  »  1,  Equation  (3)  describes  an  exponential 
dependence  of  velocity  on  the  force: 


In  a  future  publication,  the  exponential  increase  of  the  chain  slip 
rate  through  entanglements  with  increasing  tensile  force  will  be 
related  to  the  shape  of  a-e  curves  observed  for  polymer  materials. 
We  speculate  that  the  yield  stress  is  related  to  crossover  of  the 
chain  slip  rate  from  Equation  (5)  to  Equation  (6). 

Following  the  above  concept  of  polymer  plastic  deformation,  we 
express  macroscopic  stress,  a,  through  microscopic  tension  force,  F, 
and  the  rate  of  macroscopic  deformation,  de/dt,  through  the 
microscopic  chain  velocity,  V,  predicted  by  Equation  (3).  Then, 
following  the  well-known  Maxwell  model  of  viscoelastic  defor¬ 
mation  with  a  specific  type  of  viscous  damper  described  by  Equa¬ 
tion.  (3),  we  obtain  the  microscopic  equation  that  relates  polymer 
deformation,  e,  with  external  stress,  <r,  is  presented  below 


da  _  de 
dt  ~  dt 


(7) 


Here  t  is  the  polymer  relaxation  time  and  E  —  da/de  is  Young’s 
modulus  for  the  polymer.  aT  =  kTIVret,  where  Vre\  is  a  typical  volume 
of  polymer  matrix  around  the  entanglement  disturbed  by  one 
elementary  act  of  monomer  transfer  through  the  entanglement. 
The  first  term  in  right  hand  part  of  Equation  (7)  is  the  rate  of  the 
stress  change  induced  by  polymer  deformation.  The  second  term  is 
rate  the  plastic  relaxation  of  the  stress.  This  term  represents  an 
extension  of  Maxwell’s  relaxation  rate,  cr/x.  In  the  case  of  low  stress, 
a  « <77,  hyperbolic  sine,  sinh,  in  the  right  hand  side  of  Equation  (7) 
is  expanded  to  reproduce  the  well-known  Maxwell  equation  for 
polymer  viscoelastic  deformation: 


^-Es]  =  JL  (8) 

The  fundamental  difference  between  Equation  (7)  and  the 
Maxwell  Equation  (8)  is  in  the  relaxation  terms  in  the  right  hand 
sides  of  Equations  (7)  and  (8).  At  large  stress  a  »  or  the  hyperbolic 
sine  asymptotically  approaches  the  exponent,  leading  to  a  sharp 
increase  in  the  stress  relaxation  rate  as  the  applied  stress  is 
increased. 


3.  Membrane  stress  model 


Equation  (3)  indicates  non-linear  dependence  of  the  chain  slip 
rate  through  the  entanglement  on  the  tensile  force  of  the  chain. 
Equation  (3)  describes  the  conventional  linear  dependence  of 
velocity  on  the  force,  when  FalkT  «  1: 


In  this  section,  we  present  the  model  that  calculates  the  stress 
induced  by  hydration/dehydration  cycling.  In  a  fuel  cell  stack,  the 
membrane  electrode  assembly  (MEA)  is  constrained  between  two 
bipolar  plates  (Fig.  3a).  The  membrane  can  neither  bend  nor  change 
length.  Fuel  is  supplied  to  the  MEA  through  fuel  gas  channels  and 
oxidant  is  supplied  through  oxidant  (air)  channels  as  indicated 
schematically  in  Fig.  3. 


Fig.  2.  (a)  Polymer  chain  segment  constrained  by  one  entanglement  and  subjected  to  the  external  force  F.  (b)  Potential  energy  of  monomer  in  entanglement  before  applying  the 
force  (in  the  left  hand  side  of  the  Fig.  2b)  and  after  applying  the  force  (in  the  right  hand  side  of  the  Fig.  2b). 
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Fig.  3.  Schematic  figure  of  membrane  in  fuel  cell  (counter-flow  configuration)  con¬ 
strained  between  two  bi-polar  plates  (a)  and  crack  formation  (b).  The  wet  fuel  exits 
from  the  anode  gas  channel  (top  channel).  Initially,  wet  air  is  fed  into  cathode  gas 
channel  (bottom  channel).  When  cathode  side  feed  gas  (air)  is  switched  from  wet  to 
dry,  cathode  side  of  the  membrane  dehydrates  and  attempts  to  shrink.  The  anode  side 
of  the  membrane  remains  wet  and  swollen. 


To  illustrate  the  mechanism  of  membrane  failure  due  to  RH 
cycling  in  the  gas  channels  we  consider  the  following  experiment. 
Initially,  wet  fuel  is  fed  into  the  anode  gas  channel  (top  channel 
shown  in  Fig.  3a)  and  wet  air  is  fed  into  the  cathode  gas  channel 
(bottom  channel  shown  in  Fig.  3a).  Both  sides  of  the  membrane  are 
equilibrated  with  wet  gas  until  the  equilibrium  water  content  Ao  is 
reached  and  all  mechanical  stresses  in  the  membrane  relax  to  zero. 
The  cathode  side  of  the  membrane  dehydrates  and  attempts  to 
shrink  when  feed  gas  (air)  is  switched  from  wet  to  dry  in  the 
cathode  gas  channel.  Mechanical  constraints  which  prevent  the 
membrane  from  shrinking,  cause  in-plane  tensile  stress,  Oyy,  in  the 
membrane.  Tensile  stress  results  in  a  crack  formation  on  the  RFI 
cycled  side  of  the  membrane  (Fig.  3b).  Further  propagation  of  the 
crack  in  the  through-plane  direction  causes  membrane  mechanical 
failure.  These  reasoning  show  that  the  cracks  first  appear  in  the 
parts  of  the  membrane  subjected  to  the  highest  RH  cycling 
amplitude,  that  is  near  the  inlet  of  the  cathode  gas  channel,  in  the 
configuration  studied  in  this  article.  Therefore,  one  should  pay 
special  attention  to  this  region  of  the  membrane. 

There  are  many  models  in  the  literature  that  describe  water 
transport  in  the  fuel  cell  [12—14],  In  general,  the  main  problem  of 
water  transport  in  fuel  cell  is  caused  by  possible  presence  of  liquid 
water  in  the  gas  channels  and/or  in  the  GDL.  For  example,  water 
accumulation  and  removal  from  the  GDL  can  have  a  long  transient 
time  [39]  and  substantially  complicate  the  water  transport 
modeling.  Fortunately,  we  are  interested  in  the  region  close  to  the 
air  inlet,  where  there  is  no  liquid  water  in  the  GDL.  Therefore,  we 
can  use  simple  water  transport  model,  where  only  water  vapor  is 
predominant.  Certainly,  this  simple  model  can  lead  to  substantial 
inaccuracy  in  water  distribution  far  from  the  air  inlet,  but  this 
possible  inaccuracy  has  negligible  impact  on  the  membrane 
degradation  rate  near  the  gas  channel  inlet.  In  the  anode  gas 
channel  we  always  assume  wet  fuel  at  the  inlet,  so  that  the 
membrane  interface  on  the  anode  side  is  assumed  to  be  at 
RH  =  100%  and  for  simplicity,  we  did  not  solve  water  transport 
equation  on  the  anode  side. 


Within  the  cathode  side  of  the  membrane,  the  most  dry  zone  is 
the  membrane  interface  over  the  gas  channel  (near  the  inlet).  This 
region  is  anticipated  to  have  the  most  severe  degradation  condi¬ 
tions  and  therefore,  is  the  region  of  focus  for  this  article.  Further¬ 
more,  the  region  of  GDL  over  the  channel  rib  is  anticipated  to  have 
better  hydration  and  hence  is  not  considered  in  this  study.  In  this 
gas  channel  part  of  the  fuel  cell  we  assume  one-dimensional  water 
transport  through  the  GDL.  2  s-1,  it  by  many  orders  of  magnitude 
less  that  the  cycling  period.  Therefore,  for  our  calculations  we  took 
into  account  the  diffusion  resistance  of  the  GDL  in  steady  state 
approximation,  without  solving  transient  equations  for  water  vapor 
transport  in  GDL.  For  water  distribution  in  the  cathode  gas  channel, 
RH(y,t),  we  used  standard  transport  equations  [14], 

To  calculate  water  distribution  in  the  membrane  cross-section  at 
a  fixed  coordinate  y  along  the  air  channel  we  use  RH(y,t),  the 
relative  humidity  distribution  function.  Here  we  assume  that  the 
equilibration  time  of  water  in  the  membrane  is  much  smaller  than 
RH  cycling  period  in  the  gas  channels.  Simple  estimate  of  the 
diffusion  time  for  water  transport  through  the  GDL  is  0.5  ms 
( i~L^DL/Dgas  shows  that  for  typical  GDL  thickness  100  pm  and  gas 
diffusivity  Dgas  ~  0.2  cm2  s_1 ).  The  cycling  time  for  the  reactant  flow 
is  assumed  to  be  in  seconds  or  minutes.  Therefore,  the  water  at  the 
membrane  interfaces  is  in  thermodynamic  equilibrium  with  vapor 
in  the  gas  channels: 

X(x  =  0,y,t)  =  Xeq(RHc(y,t)) 

A(x  =  Lm,y,t)  =  Xeq(RHA(y,t)) 

Here  x  is  the  through-membrane  coordinate,  Aeq(RH)  is  the 
equilibrium  water  content  in  the  membrane  (an  experimentally 
measured  function  of  RH),  RHc  is  the  gas  RH  in  the  cathode  gas 
channel,  and  RH&  is  the  gas  RH  in  the  anode  gas  channel.  In  the 
current  paper,  we  utilize  the  following  conventional  diffusion 
equation,  with  a  water  content  dependent  diffusion  coefficient,  to 
calculate  water  profile  in  the  membrane: 

ffi-s(D«5)=°  (,0> 

Here  D(X)  is  the  experimentally  measured  water  diffusion 
coefficient  in  the  membrane  (see  Fig.  5).  Initial  conditions  for 
Equation  (10)  can  be  chosen  arbitrarily,  because  after  several  cycles 
the  membrane  evolves  to  a  new  quasi-equilibrium  state  governed 
by  cyclic  conditions,  and  the  initial  state  becomes  irrelevant.  We 
chose  the  initial  conditions  A(x,0)  =  Ao.  Solution  of  Equation  (10) 
with  the  boundary  conditions  Equation  (9)  gives  the  water  distri¬ 
bution  in  the  membrane,  X(x,y,t).  Equation  (10)  was  solved 
numerically  using  the  standard  MATLAB  tool. 

Below,  we  calculate  the  membrane  stress  at  the  membrane 
boundary  adjacent  to  the  air  inlet  (fixed  coordinate  y)  where  RH 
cycling  and  stress  amplitudes  are  maximal  and  limit  membrane 
lifetime.  We  start  with  a  calculation  of  the  linear  elastic  response  of 
the  membrane  to  small  changes  in  A  and  subsequently  model  the 
more  general  viscoelastic  case. 

The  macroscopic  state  of  the  membrane  is  determined  by  two 
parameters:  the  membrane  length,  Lo,  and  water  content,  Ao.  We 
assume  that  at  equilibrium,  the  stress  in  the  membrane  is  equal  to 
zero.  The  stress  in  the  membrane  can  be  generated  by  deviation  of 
the  membrane  length  from  equilibrium,  A L  =  L  -  L0,  at  constant 
water  content.  Also,  the  stress  can  be  generated  by  deviation  of 
water  content  from  equilibrium,  A A(x,t)  =  A(x,t)  -  Ao,  in  the 
geometrically  constrained  membrane  (at  constant  length).  By 
analogy  with  linear  elasticity  theory  with  thermal  expansion,  we 
calculate  the  linear  response  of  the  membrane  stress  to  small 
deviation  from  equilibrium  as  follows: 
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Stress  (MPa) 

Fig.  4.  Fitting  Tang’s  experimental  data  [3]  (triangles)  by  Equation  (22)  (solid  line). 
Fitting  parameters  values  are:  oT  =  1  MPa  and  e0  =  3.7  x  10-6. 


oyv  =  E(eyy-aM(x,t))  (11) 

Here  eyy  —  AI/Io  and  E  is  the  Young  modulus  of  the  membrane. 
Swelling  coefficient  a  determines  the  length  change,  Ay,  of 
unconstrained  (zero  a)  membrane  caused  by  the  change  of  the 
membrane  water  content  AA:  A y/y  =  aAA.  In  ID  approximation,  the 
local  stress  in  the  membrane  depends  on  the  change  of  the  local 
membrane  hydration,  (A (x,t)  -  Ao). 

In  the  membrane,  stress  can  be  induced  by  changes  in  defor¬ 
mation,  water  content  and  simultaneous  plastic  relaxation.  The 
total  rate  of  the  stress  change  is  the  sum  of  the  stress  change  rate 
caused  by  membrane  deformation  and  the  stress  change  rate 
caused  by  the  change  of  membrane  water  content 

dffyy  _  SOyyl  Qgyy  I 

~dtT  ~  dt  L  ^aT|A  K  ’ 

According  to  Equation  (11 )  the  rate  of  the  stress  change  induced 
by  the  change  of  the  rate  of  deviation  of  water  content  in  con¬ 
strained  system  with  constant  deformation  and  without  plastic 
relaxation  is 


3<7yy| 

Drjk 


(12a) 


<,2b» 

Substituting  Equations  (12a)  and  (12b)  into  the  right  hand  side 
of  Equation  (12)  we  obtain  the  equation  for  the  total  rate  of  stress 
change  in  the  membrane 

dayy(x,t)  =  £deyy(x,t)\  +Ead\X(x,t)_aTs.nh(oyy(x,t)\ 

dt  St  dt  x  V  <77  J  K  ' 

The  first  term  in  the  right  hand  side  of  Equation  (13)  is  equal  to 
zero  in  a  constrained  membrane.  Local  water  content,  A(x,t),  is 
calculated  by  Equation  (10)  with  boundary  conditions  Equation  (9). 
Equation  (13)  relates  the  stress  in  the  membrane  with  water 
content  in  the  membrane  for  arbitrary  RH  cycling  protocol.  In 
a  working  fuel  cell,  the  membrane  is  typically  compressed  more 
under  the  ribs  of  the  bipolar  plates,  resulting  in  additional  stress  in 
the  membrane  polymer  under  the  rib  area.  However,  based  on  the 
modeling  insights,  this  additional  stress  is  anticipated  only  for 
a  limited  number  of  early  hydration  cycles.  With  every  hydration 
cycle  the  plastic  flow  of  membrane  material  would  reduce  addi¬ 
tional  stress.  That  probably  can  be  modeled  within  the  framework 
that  is  presented  in  the  manuscript,  but  would  also  require  some 
modifications. 

Equation  (13)  was  solved  numerically  with  parameters  specified 
in  the  Table  1.  The  calculated  stress  as  a  function  of  time  is  pre¬ 
sented  in  Fig.  7.  The  results  indicate  that  the  membrane  stress 
becomes  a  periodic  function  of  time  after  several  cycles  if  RH  is 
a  periodic  function  of  time.  To  derive  an  analytical  equation  for  the 
periodic  stress  we  utilize  the  separation  of  slow  variable  method. 
The  plastic  (irreversible)  deformation  of  the  membrane  during  one 
cycle  is  small  because  the  cycling  period  is  much  smaller  than  the 
membrane  relaxation  time,  Tq,c  <<  x.  Therefore,  we  can  consider 
the  membrane  as  an  elastic  media  during  one  cycle.  However,  the 
membrane  slowly  approaches  a  new  equilibrium  state  during  each 
cycle.  After  a  long  time,  t »  i,  the  membrane  reaches  the  new  state 
with  new  equilibrium  water  content,  A(x).  During  the  cycle,  A 
oscillates  below  and  above  A(x),  At  A  =  A(x)  the  membrane  is  not 
stressed.  The  in-plane  stress  in  the  periodical  regime  is  driven  by 
the  deviation  of  A  from  A(x)  and  is  calculated  by  following  equation: 


The  second  term  in  Equation  (12)  is  presented  by  Equation  (7)  =  qa(X(X  q  _  J(xyj 


(14) 


Fig.  5.  Dependence  of  diffusion  coefficient  in  Nation  membrane  on  water  content  A 
(UTRC  internal  experimental  data).  There  are  not  available  data  in  range  from  A  =  2.5  to 
A  =  15  so  linear  approximation  is  used  in  this  region. 


Table  1 

Model  parameters  and  experimental  conditions 
Parameter  Description 

Operational  conditions 
TcyC  RH  cycling  time  period 

RHmax  Maximum  relative  humidity 

in  cathode  channel 

RHmin  Minimum  relative  humidity 

in  cathode  channel 
T  Cell  temperature 

Membrane  properties 

E(RH,T)  Young’s  modulus  of  the  membrane 

A(a)  Membrane  water  content  as  a  function 

of  vapor  activity  [mol  H20/mol  S03] 
D(A),  Diffusion  coefficient  of  water  in  the 

membrane 

a  Dimension  change  coefficient 

oT  Fitting  parameter 

(internal  data  fit  to  Equation.  27) 

N0  Fitting  parameter 

(internal  data  fit  to  Equation.  27) 

[3  DMA  to  in-cell  calibration  parameter 


Value/Source 


30  s 
100% 


2.3 
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Fig.  6.  Model  prediction  of  dependence  of  RH  on  time  and  X  on  time  at  the  anode  side, 
at  the  middle  of  the  membrane  and  at  the  cathode  side. 


The  stress  averaged  over  time  in  the  periodic  regime  is  equal  to 
zero  for  any  x.  Averaging  Equation  (14)  over  time  and  taking 
advantage  of  the  condition  (<%(*,  t))t  =  0,  we  obtain  the  following 
expression  for  A(x): 


irreversible  elongation,  i.e.  plastic  deformation,  of  the  membrane. 
Accumulation  of  irreversible  elongation  in  the  membrane  causes 
membrane  damage  after  a  large  number  of  cycles. 

Tang  et  al.  [3]  experimentally  studied  Nation  111  membrane 
irreversible  elongation  under  cyclic  stress.  They  observed  that 
irreversible  elongation  is  accumulated  over  a  large  number  of 
cycles  and  causes  damage  of  the  membrane  even  at  relatively 
small  amplitudes  of  cyclic  stress.  They  also  demonstrated  that  at 
relatively  small  stress,  a  <  4  MPa,  the  membrane  elongation 
rate  is  a  linear  function  of  applied  stress  magnitude.  At  larger 
stress  amplitudes  of  about  a  =  6.5  MPa,  the  elongation  rate 
rapidly  increases.  We  speculate  that  such  non-linear  depen¬ 
dence  of  elongation  rate  on  stress  amplitude  is  caused  by 
exponential  dependence  of  the  stress  relaxation  rate  on  the 
stress  magnitude  predicted  by  extended  Eyring  model  presented 
in  Section  3. 

The  elongation  of  the  membrane  subjected  to  step-like  cyclic 
stress  is  a  function  of  the  stress  amplitude  a  and  the  cycle  period 
TcyC.  Applying  Equation  (7)  to  the  membrane  elongation  during 
constant  stress  hold  in  one  cycle  we  obtain: 


ei 


dpTcyc 

Be  ' 


sinhfe) 


(17) 


(E(x,t)A(x,  t))t 
A(X)  -  {E(xft))t 


(15) 


Substituting  Equation  (15)  into  Equation  (14)  we  obtain  the  final 
equation  for  the  stress  under  periodic  RH  cycling: 


Oyy(x,  t )  =  -E(x,  t)a(x(x,  t)  -  ^[^t))’^)  (16) 

Substituting  /(x,t)  calculated  from  Equation  (10)  with  boundary 
conditions  Equation  (9)  into  Equation  (16),  we  calculate  cyclic  stress 
in  the  membrane. 


Only  a  small  fraction,  y,  of  this  elongation  is  irreversible,  i.e. 
leads  to  membrane  damage.  In  the  current  work  we  assume  that  y 
depends  only  on  cyclic  period,  T^,  and  is  independent  of  a.  Also, 
we  do  not  account  for  self-acceleration  of  membrane  damage,  i.e. 
we  assume  that  irreversible  elongation  is  proportional  to  the 
number  of  cycles  N.  The  equation  for  irreversible  elongation  of  the 
membrane  per  N  cycles  is: 

e(<r,l V)  =  NeosinhQ0  (18) 

Here  eo  is: 


4.  Damage  accrual  model 

In  this  section,  we  describe  the  model  that  calculates  the 
membrane  lifetime  as  a  function  of  applied  cyclic  stress.  Using  the 
membrane  stress  model  presented  in  Section  3,  we  calculate  the 
stress  profile  in  the  membrane,  ayy(x,t),  for  a  given  water  content, 
X(x,t).  We  predict  the  membrane  lifetime  as  a  function  of 
membrane  stress  using  the  damage  accrual  model. 

The  typical  stress  in  the  membrane  under  fuel  cell  operating 
conditions  is  much  smaller  than  the  membrane  tensile  strength. 
Many  RH  cycles  are  required  to  cause  substantial  membrane 
damage.  However,  with  each  cycle,  the  cyclic  stress  causes  small 


Fig.  7.  Model  prediction  of  stress  as  a  function  of  time  under  RH  cycling  at  the  cathode 
and  at  the  anode  sides  of  membrane. 


^  =  OtT'yrtfo c)  (19) 

The  elongation  as  a  function  of  the  stress  and  the  number  of 
cycles,  e(a,N)  was  measured  experimentally  using  Dynamic 
Mechanical  Analysis  (DMA)  equipment  for  two  specific  values  of 
stress  amplitude.  Because  eo  is  proportional  to  an  unknown 
parameter  y  and  inversely  proportional  to  another  unknown 
parameter  x,  the  lumped  parameter  eo  was  used  for  fitting.  The 
model  parameters  eo  and  <tt  were  fitted  to  DMA  data.  The 
membrane  elongation  for  arbitrary  stress  amplitude  a  was  calcu¬ 
lated  from  Equation  (18). 

To  validate  Equation  (18),  elongation  per  one  stress  cycle,  ei(a), 
was  calculated  from  Tang’s  experimental  data  [3]  for  Nafion  111 
membrane  and  plotted  as  a  function  of  applied  stress,  a,  see  Fig.  4.  It 
follows  from  Equation  (18)  that 

&'s=*»sinh(£)  <2o> 

The  experimental  data  are  in  good  agreement  with  the  predic¬ 
tion  of  Equation  (20)  with  aT  =  1  MPa  and  eo  =  3.7  x  10-6  (Fig.  4). 

Membrane  failure  occurs  after  a  critical  irreversible  elongation 
(damage),  cent,  is  accumulated.  The  number  of  cycles  to  failure,  Ncrit, 
is  calculated  from  the  following  condition: 

Ncritel(a)  —  ecrit  (21) 

We  assume  that  ecnt  does  not  depend  on  stress  amplitude  and 
frequency.  Hence  can  be  measured  at  a  specific  stress  in  a  DMA  test. 
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The  number  of  cycles  to  failure,  Nqma. 
calculated  by  equation 


sinh(<r/oy) 


out-off-cell  tests  is 


(22) 


N0m 


eo 


^rTcycl  (Jcyc) 


(23) 


We  use  Nq  as  a  lumped  fitting  parameter  obtained  from  a  DMA 
test.  To  obtain  the  values  of  No  and  or  we  performed  stress  cycling 
for  two  stress  amplitudes,  a,  and  fitted  the  experimental  number  of 
cycles  to  failure,  JVOTt,  by  Equation  (22).  Equation  (22)  predicts 
membrane  lifetime  in  a  DMA  test  with  arbitrary  stress  amplitude  a 
using  the  obtained  values  of  No  and  oj  from  experimental  data. 
Analytical  dependence  of  these  parameters  on  temperature  and  RH 
will  be  published  later. 

Fuel  cell  conditions  under  RH  cycling  differ  from  DMA  test 
conditions  because  in  DMA  testing  the  stress  is  uniformly  distrib¬ 
uted  through  the  membrane  cross-section,  while  in  the  cell,  the 
stress  in  the  membrane  changes  substantially  from  the  anode  to  the 
cathode  side.  We  incorporated  calibration  parameter  p  into  the 
current  version  of  the  model  in  order  to  use  Equation  (22)  for  in-cell 
membrane  lifetime  prediction.  This  parameter  is  used  for 
renormalization  of  stress.  The  membrane  damage  in  DMA  testing 
with  stress  amplitude  a  is  equivalent  to  membrane  damage  in  a  fuel 
cell  stack  when  the  membrane  is  subjected  to  RH  cycling  that  cause 
a  stress  amplitude  of  (3  x  a. 

To  obtain  the  calibration  parameter  (1  we  performed  a  single  RH 
cycling  experiment  in  a  fuel  cell  up  to  membrane  failure  for 
particular  conditions  and  measured  the  experimental  membrane 
lifetime  Nceu.  Using  the  membrane  stress  model  we  calculated  the 
stress  amplitude  amax  for  these  particular  conditions.  Parameter  [1  is 
calculated  from  the  following  equation: 


'-■SEfcw  (24) 

Using  the  parameters  N0  and  oT  obtained  from  fitting  DMA  data 
by  Equation  (22),  and  the  parameter  |3  obtained  from  the  in-cell 
test,  we  predict  the  in-cell  membrane  lifetime  under  arbitrary  RH 
cycling  conditions. 


5.  Model  parameters 

In  this  section  we  present  input  model  parameter  values  and 
experimental  data  for  a  GORE-SELECT®  membrane  in  a  PRIMEA® 
MEA.  In  RH  cycling  experiments,  hydrogen  with  100%  RH  was  fed 
to  the  anode  channel.  Air  fed  to  the  cathode  channel  was  cycled 
from  RHmax  to  RHmin,  with  equal  time  intervals.  Experimental 
conditions  are  summarized  in  Table  1.  Mechanical  properties  of 
GORE-SELECT  membranes  were  obtained  from  available  literature 
and  summarized  in  Table  1.  DMA  tests  required  for  model  cali¬ 
bration  were  also  performed,  using  Q800  DMA,  from  Thermal 
Analysis  Instruments. 


6.  Model  results 

In  this  section  we  report  the  result  of  modeling  of  GORE-SELECT 
membrane  damage  and  lifetime  under  in-cell  RH  cycling  condi¬ 
tions.  The  model  components,  i.e.  fuel  cell  model,  membrane  stress 
model  and  stress  accrual  model  were  incorporated  into  the  MAT- 
LAB  code.  The  model  parameters  and  test  conditions  are  summa¬ 
rized  in  Section  5. 


The  relative  humidity  variation  in  the  cathode  gas  channel 
during  the  cycle  was  calculated  from  the  fuel  cell  model  and  used  as 
the  input  for  the  membrane  stress  model.  The  calculated  relative 
humidity  variation  is  shown  in  Fig.  6.  The  function  with  fast 
exponential  relaxation  was  used  instead  of  step-wise  function  for 
better  convergence  of  the  numerical  solution.  The  calculated 
dependence  of  A  on  time  at  three  points  across  the  membrane  (at 
the  cathode  side,  at  the  middle  of  the  membrane  and  at  the  anode 
side)  is  also  shown  in  Fig.  6.  This  dependence  is  calculated  from 
water  diffusion  equation  in  the  membrane  Equation  (10)  with 
boundary  conditions  Equation  (9)  and  initial  condition  A(x,0)  =  11. 
Fig.  9  indicates  that  the  membrane  water  content  A(t)  follows  the 
humidity  cycle  RH(t).  There  is  no  time  lag  because  water  diffusion 
time  in  the  membrane  tdif  -  h2/D  ~  0.2  s  is  much  smaller  than 
cycling  period  TcyC  —  30  s.  Here  h  —  18  pm  is  membrane  thickness 
and  D  =  0.5  x  10-6  cm2  s-1  is  the  average  water  diffusion  coefficient 
in  the  membrane. 

Water  content  in  the  membrane,  A(x,t),  is  utilized  in  the 
membrane  stress  model  to  predict  the  stress  a(x,t)  distribution  in 
the  membrane.  The  model  prediction  of  the  stress  in  the 
membrane  by  Equation  (13)  as  a  function  of  time  is  shown  in 
Fig.  7  at  the  cathode  and  anode  sides  of  the  membrane.  The 
membrane  reaches  quasi-steady  state  regime  after  t  »  x,  and 
stress  becomes  a  periodic  function  of  time,  oscillating  around  the 
equilibrium  value.  The  model  prediction  for  the  stress  in  the 
periodic  regime  as  a  function  of  time  calculated  by  Equation  (16) 
is  shown  in  Fig.  8.  Maximal  stress  is  imposed  in  the  region  with 
maximal  humidity  variation,  i.e.  at  the  cathode  side  of  the 
membrane.  At  the  anode  side  of  the  membrane  the  stress  is  equal 
to  zero  because  there  is  no  humidity  cycling.  One  can  see  that  the 
stress  calculated  by  Equation  (13)  at  large  time  is  in  good  agree¬ 
ment  with  the  periodic  stress  shown  in  Fig.  8,  calculated  analyt¬ 
ically  from  Equation  (16). 

Figs.  7  and  8  indicate  that  after  time  t  larger  than  relaxation  time 
t  the  stress  at  the  cathode  sides  oscillates  around  the  equilibrium 
stress  (Teq  —  0.  The  membrane  at  the  cathode  side  is  subjected  to 
tensile  stress  when  2(0,  t)  <  A(x  =  0),  and  to  compressive  stress 
when  2(0,  tr)  >  A(x  =  0). 

The  maximal  value  of  the  stress  calculated  by  the  membrane 
stress  model  is  used  in  the  damage  accrual  model  to  predict 
membrane  lifetime  by  Equation  (24).  As  discussed  in  Section  6,  this 
equation  contains  three  unknown  parameters,  a-p,  N0  and  p. 
Parameters  cjt  and  No  are  calculated  by  fitting  DMA  data  (on  the 
number  of  cycles  to  failure  at  a  given  stress  level)  by  Equation  (22). 
Calibration  parameter  P  is  calculated  from  an  in-cell  lifetime 
experiment  under  specific  RH  cycling  conditions.  The  in-cell 
membrane  stress  was  calculated  by  the  membrane  stress  model. 


Fig.  8.  Model  prediction  of  stress  in  quasi-steady-state  regime  at  the  cathode  side,  at 
the  middle  of  the  membrane  and  at  the  anode  side. 
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Life-time,  cycles 


RH  min 

Fig.  9.  Calculated  membrane  lifetime  (the  number  of  cycles  to  failure)  as  a  function  of 
minimal  RH  at  the  cathode  side  of  membrane.  Anode  side  of  membrane  maintains 
100%  RH. 


Using  the  model  parameters  calculated  from  DMA  testing  and 
in-cell  RH  cycling,  we  predicted  GORE-SELECT  membrane  lifetime 
under  arbitrary  RH  cycling  conditions.  Predicted  membrane  life¬ 
time  (number  of  cycles  to  failure)  as  a  function  of  minimal  RH  in  the 
cathode  gas  channel,  RHmjn,  is  plotted  in  Fig.  9.  The  anode  gas 
humidity  was  assumed  to  be  100%.  The  maximal  RH  in  the  cathode 
gas  channel,  RHmax,  of  100%  was  also  assumed.  At  high  RH  cycling 
amplitude  (low  RHmin)  the  model  predicts  exponential  decrease  of 
membrane  lifetime  with  increase  of  RH  cycling  amplitude.  Stress 
amplitude  is  approximately  proportional  to  RH  cycling  amplitude, 
and  the  number  of  cycles  to  failure  depends  exponentially  on  the 
stress  for  large  stress  according  to  Equation  (22).  This  results  in 
exponential  dependence  of  the  membrane  lifetime  on  RH  cycling 
amplitude  for  large  values  of  amplitude  (small  value  of  RHmin). 

7.  Conclusions 

A  modeling  framework  was  developed  that  predicts  the  lifetime 
of  PEM  fuel  cell  membranes  subjected  to  hydration  cycling.  The 
developed  model  predicts  membrane  lifetime  as  a  function  of  RH 
cycling  amplitude  and  membrane  mechanical  properties. 
Membrane  failure  in  the  fuel  cell  is  caused  by  damage  accumulation 
under  cyclic  stress  in  the  membrane  subjected  to  cyclic  hydration/ 
dehydration.  A  fuel  cell  membrane  is  typically  subjected  to  hydra¬ 
tion/dehydration  under  fuel  cell  conditions.  One  side  of  the 
membrane  dehydrates  and  attempts  to  shrink  when  the  feed  gas 
(air)  is  switched  from  wet  to  dry  in  the  cathode  gas  channel. 
Mechanical  constraints  imposed  by  bipolar  plates  prevent  the 
membrane  from  shrinking.  This  causes  in-plane  tensile  stress  in  the 
membrane.  Tensile  stress  causes  a  crack  formation  at  RH  cycled  side 
of  the  membrane  (Fig.  3).  Further  propagation  of  the  crack  in  the 
through-plane  direction  causes  membrane  mechanical  failure. 

The  modeling  framework  consists  of  three  components: 
a  model  of  RH  distribution  in  gas  channels,  a  membrane  stress 
model,  and  a  damage  accrual  model.  Several  models  of  RH  distri¬ 
bution  in  gas  channels  are  available  in  literature  [8—12]  and  we  do 
not  discuss  this  model  component  in  the  current  paper.  In  the 
current  version  of  the  model,  we  assume  equilibrium  of  water  in 
the  membrane  with  the  vapor  at  the  membrane/gas  interface  and 
use  the  conventional  diffusion  equation  for  calculation  of  water 
content  in  the  membrane  bulk,  with  the  water  diffusion  coefficient 


dependent  on  water  content,  D(X ).  The  membrane  stress  model 
calculates  the  stress  in  the  membrane  caused  by  membrane  cyclic 
swelling/shrinking  under  RH  cycling  conditions.  The  local  stress  in 
the  membrane  is  caused  by  the  change  in  the  local  hydration  level. 
The  damage  accrual  model  predicts  the  number  of  cycles  to  failure 
for  the  membrane  under  applied  cyclic  stress.  The  input  for  the 
damage  accrual  model  is  a  maximal  stress  in  the  membrane 
calculated  by  the  membrane  stress  model  and  experimental 
membrane  lifetimes  in  DMA  tests  for  two  cycling  amplitudes.  The 
current  version  of  the  model  also  contains  one  calibration  param¬ 
eter  obtained  from  an  in-cell  RH  cycling  experiment  with  specific 
RH  cycling  conditions. 

The  model  was  utilized  for  in-cell  lifetime  predictions  of  GORE- 
SELECT  membranes.  The  membrane  mechanical  properties  and 
swelling  coefficient  were  obtained  from  the  literature.  Controlled 
tests  such  as  DMA  and  in-cell  RH  cycling  were  carried  out  in-house 
to  fit  few  model  parameters.  After  calibration,  the  model  predicts 
membrane  lifetime  in  fuel  cells  under  arbitrary  RH  cycling  condi¬ 
tions.  The  calculated  membrane  lifetime  (number  of  cycles  to 
failure)  as  a  function  of  minimal  RH  in  cathode  gas  channel,  RHmjn, 
is  plotted  in  Fig.  9.  At  high  RH  cycling  amplitude  (low  RHmin)  the 
model  predicts  exponential  decrease  of  membrane  lifetime  with 
increase  of  RH  cycling  amplitude.  Stress  amplitude  is  approxi¬ 
mately  proportional  to  RH  cycling  amplitude,  and  the  number  of 
cycles  to  failure  depends  exponentially  on  stress  for  large  stress 
magnitude  according  to  Equation  (22).  This  results  in  exponential 
dependence  of  membrane  lifetime  on  RH  cycling  amplitude  for 
large  values  of  amplitude  (small  value  of  RHmin). 
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